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ABSTRACT 

We present a set of 3-dimensional, radiation-magnetohydrodynamic calculations of the gravitational 
collapse of massive (300 M©), star- forming molecular cloud cores. We show that the combined effects of 
magnetic fields and radiative feedback strongly suppress core fragmentation, leading to the production 
of single star systems rather than small clusters. We find that the two processes are efficient at 
suppressing fragmentation in different regimes, with the feedback most effective in the dense, central 
region and the magnetic field most effective in more diffuse, outer regions. Thus, the combination of 
the two is much more effective at suppressing fragmentation than either one considered in isolation. 
Our work suggests that typical massive cores, which have mass-to-flux ratios of about 2 relative to 
critical, likely form a single star system, but that cores with weaker fields may form a small star cluster. 
This result helps us understand why the observed relationship between the core mass function and 
the stellar initial mass function holds even for ~ 100 Mq cores with many thermal Jeans masses of 
material. We also demonstrate that a ~ 40 AU Keplerian disk is able to form in our simulations, 
despite the braking effect caused by the strong magnetic field. 

Subject headings: ISM: clouds — radiative transfer — stars: formation — stars: mass function — 
turbulence — (magnetohydrodynamics:) MHD 



1. INTRODUCTION 

Massive stars, which have mass > SMq, make up < 1% 
of the total stellar population, but their numbers belie 
their impact. Both the total luminosity and the ioniz- 
ing luminosity of a star are highly super-linear functions 
of mass. Thus, massive stars have a much stronger im- 
pact on their birth environments than low-mass stars do. 
Since most stars form in clusters that contain at least one 
early O star, massive stars have an important impact on 
the formation of their low-mass neighbors, whether by 
altering the thermal properties of their parent clumps by 
heating the dust, or by destroying them outright via pho- 
toionization. The latter process is so bright that it allows 
observation of the star formation rate in other galaxies. 
Finally, massive stars end their lives in supernova ex- 
plosions, which produce heavy elements and add large 
amounts of energy to the interstellar medium (ISM), con- 
tributing to the driving of its turbulence on large scales. 
Understanding the life cycle of massive stars from their 
births to their deaths is thus an important problem for 
many branches of astrophysics. 

Unfortunately, the first stage of this process - the birth 
of massive stars - remains an incompletely understood 
problem. Observationally, regions of massive star for- 
mation in our own galaxy tend to lie farther away from 
Earth than regions of low-mass star formation, meaning 
that observers have not yet been able to probe the for- 
mation process for high-mass stars at the same level of 
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detail as they have for low-mass stars. Theoretically, the 
central difficulty is the large number of mutually inter- 
acting physical processes involved. Massive stars form 
out of a supersonically turbulent, self- gravitating fluid 
with dynamically significant magnetic fields. Massive 
protostars also deeply impact their surroundings as they 
form through a variety of feedback processes, including 
magnetically-launched outflows, radiation pressure, ra- 
diative heating, and ionization. Because of the complex- 
ity of these processes, simulations of massive star forma- 
tion are able to include at most a few of these effects 
at one time. In the past several years, there has been 
much work done on massive star formation t hat ignored 
the effects of magnetic fields, both with (e.g. j Krumholz 
et al.||2007a] |2010, 2009; Cunningh am et al.|pm| ) and 
without (e.g. Giric hidis et al.n2011( ) radiative feedback. 
There has also been much work on simulating massive 
star formation that included the magnetic fi eld, but did 
not include radiative feedback (e.g. ,Seifried et al.||2011 
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20121 [Li &: Nakamura|2006l|Wang et al.|:^(jlo||HenneCTe| 
et al.pOlip . Thus far only two published sirnulations of 
massive star formation have included both radiation and 
magnetic fields, and these provide only a limited pic- 
ture of how fragm entation in massive cores works. |Pe-| 
ters et al. (2011) treat direct stellar radiation and ion- 
ization chemistry, but neglect the dust-reprocessed ra- 
diation field, w hich is mainly responsible for regulating 
fragmentation. 'Commercon et ah] (2011) include dust- 
reprocessed light, but because they do not employ a sub- 
grid stellar model they are forced to halt their calcula- 
tions when < 1% of the core material has collapsed, and 
as a result they cannot study the fragmentation of the 
bulk of the gas. 

In this paper, we attempt to fill that gap. We present 
the results of 3-dimensional, adaptive mesh refinement 
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(AMR), radiation-magnetohydrodynamic (R-MHD) sim- 
ulations that treat the dust-processed radiation from pro- 
tostars in the flux-limited diffusion (FLD) approxima- 
tion. In particular, we focus on the fragmentation of iso- 
lated, massive cores in the relatively early stages of star 
formation - up to the point at which about 10% of the 
core gas has turned into stars. The question of how mas- 
sive cores fragment is an important one for any theory of 
star formation in which the initial mass function (IMF) 
is set in the gas phase, e.g. the turbulent fragmentation 
scenario originally laid out in Padoan & Nordlund ( 2002 ) . 
Observations of the core mass function (CMF) in galactic 
star-forming regions reveal that it looks like a sca led-up 
version of the stellar initial mass function (IMF) (|Alves 



et al. 120071 [Nutter &: Ward-Thompson|2007| [Enoch et al. 
' 2008p . Thi s relationship appears to continue even up to 
^^TI)OM0 (Reid & Wilson 2006). This correspondence 
- that the CMF has the same form as the IMF but is 
shifted up in mass by a factor of ~ 3 - has a natural 
explanation if massive cores do not fragment strongly as 
they collapse, but instead simply convert ~ 1/3 of their 
mass into single massive stars or systems. 

The purpose of this paper is to address the question 
of how massive cores fragment via direct numerical sim- 
ulation. Our outline is as follows: in section 2, we de- 
scribe our numerical setup, including the equations and 
algorithms used as well as our initial and boundary condi- 
tions. In section 3, we present our results, focusing on the 
evolution of our cores over a period of 0.6 mean-density 
free-fall times. In section 4, we discuss our results, in 
which the magnetic field and the radiative transfer to- 
gether have a significant impact on the fragmentation 
of the cores in a way one would not predict from either 
process considered in isolation. We summarize our con- 
clusions in section 5. 

2. NUMERICAL SETUP 

2.1. Equations and Algorithms 

We solve the equations of mass, momentum, and en- 
ergy conservation on a hierarchy of AMR grids. We as- 
sume that the motion of the gas is governed by the ideal 
MHD equations a nd treat the radiation us ing the mixed- 
frame approach of Krumholz et al. (2007b). At any time, 
the computational domain consists of a fluid made up 
of gas, dust, and radiation, plus some number of sink 
particles that represent stars. The fluid quantities are 
described by a vector of state variables (p, pv^E^B^ Eji) 
defined at every grid cell, where p is the gas density, 
pv the momentum, E the non-gravitational energy den- 
sity (i.e. the total of the kinetic, thermal, and magnetic 
energy densities), B the magnetic field, and Er the radi- 
ation energy density. The particles are characterized by 
their position a?^, momentum p^, mass M^, and luminos- 
ity Li, which is determine d via the prot ostel lar evolution 
model described in McKee & Tan (2003) and Offner et al. 



(2009). The equations governmg the evolution of the K- 
MHD fiuid-particle system are: 
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In the above equations, the total pressure Pt is Pgas 
P^/Stt, and we use an ideal equation of state, so that 
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where /c^ is the Boltzmann constant, Tg the gas tem- 
perature, p the mean molecular weight, 7 the ratio of 
specific heats, and e the thermal energy per unit mass. 
We take p = 2.33 and 7 = 5/3, appropriate for molec- 
ular gas of solar composition that is too cold to store 
energy in rotational degrees of freedom. The corre- 
sponding value for the gas's specific heat capacity is 
Cv = ks/i^y - l)pmB_ ~ 5.3 X lO'^ erg g"^ K"^ 

The summations in the gas-sink interaction terms are 
taken over all the particles in the domain, and W{x — Xi) 
is a weighting kernel that distributes the transfer of mass, 
momentum, and energy over a radius of 4 fine-level cells 
around sink particle i. The values for M^, p^, and Si, 
or the rates of mass, momentum, and energy transfer 
between the sink particles and the fluid, are computed by 
fitting the flow around each sink particle to a magnetized 
Bondi-Hoyle flow; see Lee et al. (2013, in preparation) for 
details. The star particle states themselves are updated 
according to the following equations: 



dt 
d 



Mi = Mi. 
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Pi = -MiV(t)^pi. 
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Because our sink particle algorithm destroys information 
about the fluid flow inside the 4 fine cell accretion zone 
around each particle, we are not able to properly fol- 
low the dynamics of particles that pass within that dis- 
tance of each other. We therefore adopt the following 
criterion to handle mergers between sink particles that 
pass within one accretion radius of each other (40 AU in 
most of the simulations presented here): we merge the 
two sinks together only if the smaller sink is less than 
O.O5M0 in mass. This threshold roughly corresponds 
to the mass at which second collapse occurs (Masunaga| 
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et al.||1998l [Masunaga fc Inutsuka||2QQQD . Before that 
pomt, sink particles represent hydrostatic cores of several 
AU in size, which could be expected to merge together. 
After that point, they have collapsed down to roughly 
solar size scales, and will not necessarily merge simply 
because they pass within 40 AU of each other. 

The gravitational potential (j) in the above expressions 
obeys the Poisson equation with a right-hand side that 
includes contributions from both the fluid and the star 
particles: 



p^^Mi5{x-Xi) 
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where G is the gravitational constant. 

The radiation-specific quantities are the speed of light 
c, the comoving frame specific Planck- and Rosseland- 
mean opacities kqr and /^op, and the Planck function 
Bt = canTg /{An), where an is the radiation constant. 
Finally, the flux limiter A and Eddington factor R2 are 
two quantities that enter the flux-limited diffusion ap- 
proximation we use to c ompute the radiative transfer. In 
this work, we adopt the Levermore & Pomraning (1981) 
approximation: 
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We obtain the dust opacities a>:op a nd /^qr from a 
iecewise-linear fit to the models of 'Semenov et al. 



2003); see Cunningham et al. (2011) for the exact func 



tional form. 

We solve the above equations using a new version of our 
astrophysical AMR code Orion, which allows us to si- 
multaneously include the magnetic field and the radiative 
feedback. Orion solves the above equations in a number 
of steps, which we summarize below. First, we solve the 
ideal MHD equations by themselves (Equation (|4|, the 
first two terms of Equation ([T]) and the first three terms 
of Equations ^ and ([3|) usmg a Godunov- type scheme 



with the HLLD approximate Riemann solver (Miyoshi 
& Kusano||2005|). Specifically, we use the dimensionally 



Li et al. (2012), which makes use of the uni- 



unsplit, A MK Constraine d Transport (CT) scheme de- 
scribed in^ 

grid CT sch eme Irom the open-so urce astrophysical MHD 
code Pluto (Mignone et al.|2012). This portion of the up- 



date algorithm uses a lace-centered representation for the 
magnetic field and we use the Chombo AMR library 
to provide support for the face-centered fields. Ne xt, we 
i ncorp orate self-g r avity in the manner of T ruelove et al.| 
(1998) and Klein (1999). To solve the Poisson equation 
(Equation (|lo|)J, we use an iterative multigrid scheme 
also provide3~Dy Chombo. In the third step, we update 
Equations ([2|, ([3|, and ([5| for the radiative terms us- 
ing the opera tor-split approach described in Krumholz| 
et al. ( 2007b| ). Briefly, this technique first solves the ra- 
diation pressure, work, and advection terms explicitly, 
and then implicitly updates the gas and radiation energy 
densities for the terms involving diffusion and the emis- 
sion/absorption of radiation. This update is handled by 



the iterative process described in Shestakov et al. (2005), 
which uses psuedo-transient continuation to reduce the 
number of iterations required for convergence. We then 
complete the update cycle by calculating the new sink 
particle states using the above equations and comput- 
ing their interactions with the fluid using the algorithms 
described in Lee et al. (2013, in preparation). 

Finally, we point out some important numerical 
caveats: our treatment of the radiation in this work fo- 
cuses on the diffuse, dust-processed component of the 
radiation field, and it treats that radiation as gray. Mas- 
sive stars, however, put out large numbers of ionizing 
photons, and these photons have a dramatic impact on 
the surrounding environment. Furthermore, treating the 
diffuse component of the field as gray and ignoring the 
direct component of the non-ionizing radiation both lead 
us to underestimate the radiation pressure force by a fac- 
tor of a few (Kuiper et al.„2011) . However, since both 
of these effects are most significant for stars more mas- 
sive than ~ 2OM0, and since our conclusions are mainly 
based on the evolution of the cores prior to the most mas- 
sive star reaching that point, we do not believe that our 
qualitative conclusions will be significantly altered by a 
more accurate treatment of the radiative transfer. We 
have also not included the effects of protostellar outflows 
in any of the runs in this paper. We shall do so in future 
work. 

2.2. Refinement and Sink Creation 

The computational domain is a cube with side Lbox 
that is discretized into a coarse grid of Nq cells, so that 
the resolution on the coarse grid Axq = Li^q^/Nq. Our 
code operates within an AMR framework that automat- 
ically adds and removes finer grids as the simulations 
evolve. With L levels of refinement and a refinement 
ratio of 2, the resolution of the finest level is Axl is 
In this work, we have chosen these parameters 
such that Axl is 10 AU. 

Any cell that meets one or more of the following criteria 
is flagged for refinement: 

1. The density in the cell exceeds the magnetic Jeans 
density, given by 



GAxf 



0.74 



(14) 



where Cg is the isothermal sound speed, Axi the 
cell size on level /3 = Sirpcl/B'^ and Jmax is 
the maximum allowed number of magnetic Jeans 
lengths per cell, which must be small to avoid artifi- 
cial fragmentation. Throughout this work, we take 
^max = 1/8. Note that this is identical to our pre- 
vious work except for the inclusion of the magnetic 
field. Because the field provides additional support 
against collapse, we do not need to resolve the flow 
as highly in the presence of magnetic fields to pre- 
vent artificial fragmentation. For a derivation and 
numerical justification of this relation, see the Ap- 
pendix, but we note that it is roughly equivalent to 
including the magnetic energy density along with 
the thermal energy in the expression for the Jeans 
length. 

2. The cell is within 16 Axi of a sink particle. 
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3. The gradient in the radiation energy density ex- 
ceeds 

\/ER>0.2b^. (15) 

^ Axi ^ ^ 

This procedure is repeated recursively until the final level 
is reached. At that point, if there are still any cells on the 
finest level that exceed the magnetic Jeans density, then 
the excess matter is removed from the cell and placed 
into a new sink particle, which then evolves according to 
the algorithm in section 2.1 above. Taken together, these 
three conditions ensure that the regions where star for- 
mation is happening are always tracked with the highest 
available numerical resolution. 

The application of these criteria to simulations of self- 
gravitating, isothermal gas requires special care, because 
such simulations have a fundamental problem: They do 
not converge. Isothermal gas tends to produ ce long, thin 
filaments, which do not fragment st rongly (Inutsuka & 
Miyama 1992 Truelove et al.|p!998 ) and are thus non- 



trivial to decompose mto pomt particles. Convergence 
studies by Boss et al. (2000) and Martel et al. (2006) 
suggest that there is no well-defined, converged solution 
for fragmentation and sink particle creation in this case, 
because the correct solution is collapse to singular fila- 
ments rather than singular points. As a result, for any 
choice of the finest resolution, application of the Truelove 
criterion to a collapsing isothermal gas will result in pro- 
ducing artificial fragments at the finest grid scale. This 
does not mean that all fragmentation in isothermal sim- 
ulations is artificial: As we shall see below, our isother- 
mal simulation produces about the same total mass in 
stars and the same amount of mass in the most mas- 
sive star as our radiative simulations; on the other hand, 
it produces many more low-mass stars. In view of this 
over-fragmentation problem in isothermal simulations of 
star formation, it is essential to carry out a resolution 
study to verify that the conclusions being drawn from 
such simulations are physical and not numerical. 

Interestingly enough, while much of the fragmentation 
in isothermal simulations is ultimately caused by the nu- 
merical mesh, proper adjustment of the finest level of 
resolution may nonetheless enable isothermal simulations 
to give a qualitatively correct picture of fragmentation 
in the absence of radiative feedback. Without protostel- 
lar heating, molecular gas still becomes non-isothermal 
at some density pcrit at w hich energy can no longer be 



efficiently radiated away. |Masunaga & Inutsuka 
find that 



■ |j 

for our choice of initial temperature an( 



(11999) 

dust 

opacity, pcrit ^ lO"^'^ g cm"'^. Past that point, the ther- 
mal pressure inside the filament starts to become more 
important relative to gravity. Eventually, gravitational 
contraction begins to slow, the timescale for cylindrical 
collapse becomes large compared to that for spherical col- 
lapse, and fragmentation will occur. Unfortunately, the 
results from such a simulation cannot be validated with 
a convergence study: increasing the resolution makes the 
fragments that form smaller than appropriate for the ac- 
tual, non-isothermal case. 

We stress that the non-convergence of the number of 
fragments in isothermal simulations is not a consequence 
of our particular sink particle algorithm. Using more 
stringent sink creation criteria, like those proposed in 
[Federrath et al., (2010J , has the benefit of producing 



fewer spurious fragments, but some will still be present, 
and their properties will still ultimately be determined 
by the numerical mesh. Furthermore, one cannot get 
around this problem by suppressing sink formation en- 
tirely within filamentary structures, since once the Tru- 
elove criterion is violated the filament will fragment arti- 
ficially anyway. To get a converged answer on the number 
of fragments formed in self-gravitating, turbulent media, 
one must include some sort of deviation from isother- 
mality and a fine enough numerical mesh to resolve the 
resulting fragments. 

2.3. Initial and Boundary Conditions 

We begin with three cores that are identical except 
that we include a different combination of physical pro- 
cesses in each run. The parameters for these simulations 
are summarized in Table [TJ Run HR includes the radia- 
tive transfer physics but has no magnetic field, run BI 
has a magnetic field but no radiation, and run BR has 
both a magnetic field and the radiative transfer. For run 
BI, we have dropped Equations Q and ([5| and adopted 
the isothermal equation of state (Pgas = pcs'^) instead of 
Equation (|6|. 

With the exception of the magnetic field, our initial 
conditions a re almost identical to the those in M yers| 
et al. ( 2011|) and (with the exceptio n of the protostellar 



out flows)" Cunningham et al. (20TlJ. In all of our runs, 
we begin with an isolated sphere of gas and dust with 
mass Mc = 300 M©, radius Rc = 0.1 pc, and tempera- 
ture Tc = 20 K. The density follows a power-law profile 
proportional to r~^"^, so that the density at the edge of 
the core is 

Pedge = (16) 

The surface density of these cores, = Mc/ttR^ ~ 2.0 g 
cm~^, is chosen to resemble that observed in galactic re- 



gions o f high- mass star formation. For example, ,McKee| 
Tan (2003) inferred a mean E ^ 1 g cm~^ fro m the 



sample of high-mass clumps in Plume et al.| ( 1997). The 
corresponding mean density is p ^ 4.8 x 10~ g cm~^, or 
nn = 2.4 x 10^ H nuclei cm This value determines the 
characteristic timescale for gravitational collapse, given 

by 



Stt 
32Gp 



30.2 kyr. 



(17) 



While these initial parameters are to an extent chosen 
for computational convenience (higher densities mean 
shorter free-fall times, which mean fewer total time steps 
need to be taken) they are consistent with sub-mm in- 
terferometric observations of massive cores (Swift 2009). 
Furthermore, the r~^-^ density profile agrees with obser- 
vations of star- forming regions at the ^ 1 pc clump scale 
fBeuther e t al.| ( |2007j ), ICaselh fc Myers] ( | 1995 |, [ Mueller 
et al. (20Q 2|) and the ^ 0.1 pc core scale ( [Longmore 
|et al.||20ir "Zhang et al. 2009). Si milarly , a recent mid- 
infrared extmction study (Butle r fc Tan||20i 2) observed 
42 massive cores in 10 different IKDCs and (after enve- 
lope subtraction) reported a mean kp of ~ 1.6. They 
also report that the power-law profile was a better fit to 
their observations than the less centrally concentrated 
Bonnor-Ebert profile. 



Fragmentation in Magnetized, Massive Star-Forming Cores 



TABLE 1 
Simulation Parameters 



Name 


RT? 


M (Mo) 


R (pc) 


(Jv (km s -*■) 


(kyr) 


M/M^ 


B (mG) 




^box (pc) 




L 


AxL (AU) 


HR 


Yes 


300 


0.1 


2.3 


30.2 


oo 


0.0 


CO 


0.4 


256 


5 


10.0 


BR 


Yes 


300 


0.1 


2.3 


30.2 


2.0 


1.6 


0.05 


0.4 


256 


5 


10.0 


BI 


No 


300 


0.1 


2.3 


30.2 


2.0 


1.6 


0.05 


0.4 


256 


5 


10.0 



Note. — Col. 8: mean magnetic field in the core. Col. 9: mean plasma /5 = SirpCg/B^ in the core. Col. 10: resolution of 
the base grid. Col. 11: number of levels of refinement. Col. 12: maximum resolution at the finest level. 



Our cores are placed at the center of a cubic box with 
side length equal to 0.4 pc, so that the sides are far 
enough removed from the core that there is minimal in- 
teraction from the boundaries. The parts of the box that 
are not covered by the core are filled with a hot, diffuse 
medium with = Pedge/10 and = 200 K, so that 
the ambient medium will be in thermal pressure equilib- 
rium with the core. We set the opacity of this confining 
gas to zero so that it will not cool as the simulation pro- 
ceeds. The initial condition on Eji is given everywhere 
by aRT^^ where the radiation temperature Tr is also set 
to 20 K. 

For boundary conditions, we choose outfiow for the 
MHD update, meaning that in advancing the hyperbolic 
subsystem we set the gradients of p, pv, and B to zero 
at the domain boundary. For the radiation update, we 
use Marshak boundary conditions, meaning that the en- 
tire simulation volume is bathed in a blackbody radiative 
fiux corresponding to 20 K, while radiation generated 
within the simulation volume may escape freely. Finally, 
in solving Equation ( 10 ) for (/), we require that (/) = at 
the boundaries. 

We also give the core an initial ID velocity dispersion 
of (Jc = 2.3 km s~^, chosen to put the core into approx- 
imate virial balance. If we take the virial ratio a to be 
balRc/GMc dBertoldi fc McKee|[l992| , then a ^ 2.1. 
Thus, there is initially slightly more kinetic energy than 
gravitational potential energy in each of our cores. We 
choose a slightly super-virial value for a because we do 
not drive the turbulence by adding kinetic energy after 
the simulations begin. Although the virial parameter 
greater than unity at t = 0, it has decayed to ~ 1.0 by 
the time the simulations end. The velocities themselves 
are drawn from a Gaussian random field with power spec- 
trum P{k) cx appropriate for the highly supersonic 
turbulence found in molecular cloud cores. We include 
the perturbations in the following manner: first, we gen- 
erate a 1024^ pertur bation cube using the method of Du- 



T 



binski et al. (1995) with power on scales ranging from 



to /cmax = 512. We then place the cube over 
the simulation volume and either coarsen or interpolate 
the perturbation data so that we can represent pertur- 
bations at all levels of refinement. We have made no 
attempt to filter out compressive modes from the initial 
velocity field. The precise mixture of solenoidal and com- 
pressive components have been found to be important 
for gravitational fragment ation in unforced core collapse 
simulations (Girichidis et al. 2011} and on the overall 
rate o f star formation m simulat ions with driven turbu- 
lence (iFederrath fc Klessen|20l'2 ), but we do not explore 
this effect here. 
In our MHD runs, we also give the cores an initial mag- 



netic ffeld pointing in the z direction. The importance of 
this ffeld is best expressed in terms of the mass-to-ffux 
ratio: 

M/M^, (18) 



where 



p^ 



2^GV2 



(19) 



is the magnetic critical mass and $ is the magnetic ffux 
threading the core. Cores with > 1 are unstable 
against gravitational collapse, while cores with ja^ < 1 
are expected to be stable. Measurements of Zeeman 



splitting in both the OH molecule (Troland fc Crutcher 
"2008), which probe s densities of 1 0^ cm~'^, and the CN 



molecule (Falgarone et al. 2008), which probes higher 
densities of 10^~^ cm"'^, show that the mean value of 
jj.^ is approximatel y 2, a value su pported by theoretical 
arguments as well (McKee 1989). Note, however, that 
there may be substantial scatter in the magnetic ffeld 
strength such that many dark molecular cloud cores have 
much more supercri tical values of the mass-to-ffux ratio 
(Crutcher et al. 2010} . In this paper, we adopt ji^ = 2 for 
all of our MHD runs, and defer a more extensive param- 
eter study on the effects of the magnetic ffeld strength to 
a later work. 

In the absence of more detailed information about the 
magnetic ffeld geometry, we will assume that the spatial 
dependence of the initial B ffeld follows the cylindrically 
symmetric proffle 



B{Rz) — ^edge 



-1/2 



(20) 



where Rz is the distance to the z axis and the value of 
^edge is chosen to give the desired mean mass-to-ffux 
ratio for overall core: 



edge 



3 VGMe 

2 /u$i?2 



(21) 



For /i<|, = 2, ^edge ^ 1-2 mG. Using this form for the 
initial magnetic ffeld is clearly an idealization, but it does 
have the advantage that it 1) satisffes the condition V • 
B = 0^ and 2) ensures that the mass-to-ffux ratio in the 
central ffux tube (~ 5.6 above critical) does not greatly 
exceed the mean value for the overall core, consistent 
with the Zeeman measurements discussed above. 

Our initial conditions do not include any explicit ro- 
tation on top of the random turbulent perturbations de- 
scribed above. However, these perturbations do include 
some incidental angular momen tum. In fact, as found by 
Burkert fc Bodenheimer (2000) , Gaussian random turbu- 
lence alone may be sufficient to account for the observed 
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rotational properties of prestellar cores. When we ap- 
ply the technique in that paper to measure ^rot foi" our 
cores, we g et get Prot = 0-012, in line with the values 
observed in [Goodman et al" (|1993 ). Note, however, that 
as discussed m Dib et al. (|z010), the rotational prop- 
erties of cores nieasured in projection by observers may 
differ substantially from the actual 3D values. In fact, 
if we calculate £^rot/^grav from our initial conditions us- 
ing the full 3D velocity and density information, we get 
^ 0.002, lower than /3rot by a factor of 6. Thus, while 
the rotation in our initial conditions is consistent with 
observations, it is significantly lower than in other sim- 
ulations that impose solid-body ro tation in addition t o 
random turbulence, such as those of Seifried et al. (2012 ). 
Finally, as we do not chose the direction of the angular 
momentum vector in our cores explicitly, there was no 
imposed choice about the initial orientation of the core 
angular momentum vector L with respect to B. It turns 
out to be misaligned with the magnetic field by ^ ~ 60 
degrees. 

While the above initial conditions are clearly some- 
what artificial, they do capture the essential observed 
properties of high-mass dark-cloud cores. The most un- 
realistic aspect of our initial conditions is probably our 
imperfect treatment of the initial turbulence. While 
we include perturbations to the velocity field, there are 
no corresponding perturbations to the density at time 
t = 0. Thus, while the velocity field soon creates fila- 
mentary structures reminiscent of those expected from 
turbulence, these filaments do not have the same proper- 
ties they would in a self-consistent realization of a turbu- 



lent densit y- velocity field, as disc ussed in |Krumholz et al" 
|2012) and|Federrath & Klessen| (|2012|). [Krumholz et al. 



(2012) found that this ditt'erence can have an important 
impact on e.g. the overall star formation rate, so we 
mention it here as a caveat. Another caveat is that our 
initial velocity field does not include any infall motions 
at t = 0. This probably has the effect of encouraging 
fragmentation somewhat, since the accretion rates and 
therefore the protostellar heating rates would be higher 
if infall were included from the beginning. Ideally, one 
would generate initial conditions for massive cores from 
larger simulations at the clump scale, which would then 
contain self-consistent density perturbations and infall. 
We are considering these issues in simulations of mas- 
sive star formation at the cluster scale that are now in 
progress. The goal of this paper is to examine an ideal- 
ized case first to elucidate the underlying physics. 

We wish to emphasize that we have chosen the above 
runs to as far as possible create a controlled experiment 
where we have isolated the effect of only one physical 
process. Runs BR and BI are identical expect for the 
presence of the radiative feedback, and runs HR and BR 
are identical except for the presence of the magnetic field. 
Thus, we can isolate the effect of the radiative feedback 
by comparing the first set of runs, and the effect of the 
magnetic field by comparing the second. 

3. RESULTS 

Here, we summarize the main results of our calcula- 
tions. The simulations presented here were run on the 
NASA supercomputing platform Pleiades on 128 to 512 
processor cores and took a total of about 700,000 CPU 
hours. 



3.1. Density Structure 

The time evolution of the large-scale structure of cores 
BR, HR, and BI is shown in Figure [l] In all three runs, 
the imposed velocity perturbations create a system of 
filaments embedded within the collapsing core that feed 
gas into the central region where the massive star is form- 
ing. In the MHD runs the velocity perturbations rear- 
range the field lines so that the filaments are primarily 
perpendicular to the field. At this scale, the primary 
difference between the runs is that the filamentary struc- 
ture created by the velocity perturbations in run HR is 
much more pronounced than in either of the runs with a 
magnetic field, despite the fact that all three runs have 
the same sonic Mach number of ~ 15. There are two 
reasons for this behavior. First, even though the cores 
in runs BR and BI are highly supersonic, they are only 
marginally super- Alfvenic, with A4a ^ 1-9- The pres- 
ence of the faster magnetic signal speeds means that al- 
though shocks parallel to the magnetic field lines can be 
as strong as in run HR, flows perpendicular to the field 
that would be strong shocks in run HR are only weak 
shocks - or not shocks at all - in the other two runs. The 
overall effect is that, even ignoring gravity, the density 
contrasts imposed by the turbulence in the MHD runs 
are smaller than the hydro only run. Second, in all three 
runs, over-densities created by the turbulence can grow 
due to the self-gravity of the gas. However, in the pres- 
ence of the magnetic field, these dense regions are only 
able to grow by drawing in material along the field lines, 
whereas there is no such restriction in the hydrodynamic 
case. The combined effect is that density distribution in 
the cores at a given time is broader in run HR than in the 
other two - that is, the dense regions are more dense and 
the diffuse regions more diffuse. Finally, we note in pass- 
ing that at this scale the effect of the radiative heating 
has essentially no effect on the morphology of the core; 
the gas structure in runs BR and BI appears practically 
identical. 

The situation is different when we zoom in to show the 
central 5000 AU of the simulation volume as in Figure 
[2j where the center is defined as the location of the most 
massive star in the simulation. At this scale, we begin to 
see clear differences in the gas morphology between runs 
BR and BI. In both cases, the gas collapses into a network 
of filaments, and there is a rough correspondence between 
the filaments in BR and those in BI. However, the fila- 
ments in run BR are much fatter and more diffuse than in 
run BI. This is easily understood as a consequence of ra- 
diative heating. For an isothermal, magnetized filament 
like the ones in BI, both the magnetic and pressure forces 
scale the same way with filament size as gravity in the 
virial theorem (see the Appendix for a more detailed dis- 
cussion). Thus, either the total pressure (magnetic plus 
thermal) is initially enough to halt collapse, or else it will 
never be and the filament will collapse until something 
causes the equation of st ate to deviate from isothermality 
( [Inutsuka fc Tsuribe|200 1 ). This behavior is clearly seen 
in run BI, where the filaments contract until they reach 
the density at which our code creates sink particles. In 
run BR, on the other hand, radiative feedback from the 
central protostar has already caused the gas to become 
non-isothermal, and thus filaments close to the protostar 
stop collapsing before much sink creation takes place. 
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Fig. 1. — Column density through the simulation volume at 6 different times for runs BR (left), BI (middle), and HR (right). Projections 
are taken along the x direction, and the initial magnetic field is oriented in the positive z direction. We have set the viewing area of the 
images to be 0.3 by 0.3 pc to show the global evolution of the entire core. Star particles are portrayed as black circles, with the size of the 
circle corresponding to the mass of the star. The smallest circles represent stars with masses between O.OSM© and I.OMq. The next size 
up represents masses between I.OMq and 8.OM0, and the largest represents stars with masses greater than S.OMq. 
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Fig. 2. — Same as Figure 1, but zoomed in to show the central 5000 AU around the most massive star in each simulation. Projections 
are still taken along the x direction through the entire simulation volume. 
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Fig. 3. — Top - Face-on view of the disk in the high-resolution 
version of run BI at 0.2 %. The colors correspond to the column 
density through a sphere of radius 100 AU centered on the star 
particle. The arrows show the direction of the mean in-plane ve- 
locity of the disk gas. Bottom - the black circles show the mean 
angular velocity uj in the disk as a function of cylindrical radius 
Rz- The red line corresponds to a Keplerian profile normalized 
using the mass of the star. We have also shown the sink particle 
accretion zone in green to demarcate the radius at which our sink 
particle algorithm begins to alter the fluid properties. 

We can also isolate the effect of the magnetic field 
on the gas morphology by comparing runs BR and HR. 
There are two main differences. First, without the mag- 
netic field to help support it, the main filament of gas 
feeding the central protostar has already begun to frag- 
ment into self-gravitating, spherical "beads" by 0.3 free- 
fall times. These beads have a characteristic size of a 
few hundred AU, and are therefore well-resolved in our 
runs. The type of grid-induced filament fragmentation 
discussed in section 2.2 in the context of isothermal sim- 
ulations is thus not a concern in runs BR and HR. Sec- 
ond, beginning around the same time, we can see the 
presence of a dense, ~ 200 AU disk around the most 
massive star in run HR. This disk is centrifugally dom- 
inated with a roughly Keplerian velocity profile. We do 



not see a similar disk in either of our runs with a mag- 
netic field, at least at the 10 AU resolution of the simu- 
lations presented here. This is the well-known magnetic 
braking effect, where at = 2 the field is so efficient 
at removing angular momentum from the center of the 
core that it suppresses the formation of a Kepler i an disk 
( [ Allen et al]|2QQ3[ [Hennebelle fc Fromang||2Q08l [Mellon 
Li ' 2008[ ). However, it we repeat run Bl witn three 
more levels of refinement so that the maximum resolu- 
tion is 1.25 AU and the sink accretion radius is 5 AU, 
we do in fact begin to see a rotationally-dominated disk 
beginning around ~ 0.15 tff. By about ~ 0.2 tff, when 
the star has reached a mass of about 3.5M0, the disk 
has grown to ~ 40 AU and developed a Keplerian veloc- 
ity profile, as shown in Figure |3] This would lie entirely 
within the sink particle accretion zone in our simulation 
with 10 AU resolution, so it is not surprising that we do 
not see it there. While magnetic braking has certainly 
removed angular momentum from the material accreting 
onto the disk, allowing it to fall much closer to the cen- 
tral protostar than would be the case without a magnetic 
field, we do not find that it suppresses the formation of 
a disk entirely at high resolution. 

Several other researchers have already reported form- 
ing disks in MHD simulations of star form ation. In a 



study of magnet ic braking in low-mass cores, [Hennebelle 



fc Ciardi] ( 2009 ) found that the efficiency of magnetic 
braking depends on the angle between the initial mag- 
netic field and the core's angular momentum vector, with 
a 90 degree misalignment lowering the value of /i$ at 
which disk formation is su ppressed by a facto r of 2 — 3 
relative to the aligned case. Santos-Lima et al. ( 2012b|a ) 
studied this problem numerically as well, arguing that 
the presence of turbulence increases the rate of magnetic 
diffusion in the inertial range, allowing parcels of gas 
that have lo st magnetic ffux to fall onto a disk. Seif ried] 
et al. (2012) also found that the presence of turbulent 



perturbations reduces the efficiency of magnetic braking 
enough to form a Keplerian disk at fi = 2.6, although 
they disagree that ffux loss is involved. In our disk, we 
ffnd fi averaged over a 100 AU sphere around the most 
massive star has risen to ^ 20 by the snapshot displayed 
in Figure |3| although we have not veriffed that this is 
due to the mechanism proposed by Santos-Lima et al. 

Finally, we brieffy mention one more difference between 
our magnetic and non-magnetic runs: the presence of 
episodic outffows in runs BI and BR. Around 0.3tff, we 
begin to find material in those runs with radial veloci- 
ties of ~ 10 km s~^ away from the primary star. These 
outfiow velocities increase with time, such that by 0.6tff 
(when the primary has grown to > 2OM0) they can be 
as large as 40 km s~^, which is roughly the Keplerian 
speed at the grid scale. ~ 10 km s~^ outfiows have been 
observed previousl y in non-radiative MHD simulations of 
massive cores (e.g. Seifried et al.[[20lT Hennebelle et al. 



2011). However, because the outtlow launching mecha- 



nism is badly under-resolved in our simulations, we shall 
not discuss outfiow properties in detail here. 

3.2. Magnetic Field Structure 

Although the magnetic field lines are initially oriented 
in the z direction, this is not an equilibrium configura- 
tion, and as the simulations proceed they settle into a 
new, quasi-equilibrium "hourglass" shape shown in Fig- 
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Fig. 4. — A density slice taken through the center of the computa- 
tional domain perpendicular to the x-axis at 0.3 free-fall times. The 
y- and z- components of the magnetic field lines are over-plotted 
in white with evenly-spaced anchor points along the y-axis. 

ure [4j which rese mbles the morp hology in t he dust po- 
larization maps ofjGirart et al.|(|2009) and |Tang et al.| 
( |2QQ9| ). Here, we take a density slice through the center 
ot the domain aligned to be perpendicular to the x direc- 
tion. On top of that slice, we show the planar compo- 
nents (that is, the y- and z- components) of the magnetic 
field lines. This slice is taken from run BR at 0.3 free-fall 
times, but the overall shape of the field lines is similar at 
other times as well, provided enough time has passed for 
the initial conditions adjust to the new equilibrium. Be- 
cause the Alfven Mach number of the initial turbulence 
is ~ 2, the lines are able to be bent somewhat by the 
turbulent perturbations, but this is not a large effect. In 
the slice shown in Figure [4] we can see a dense filament 
in red, with the field lines adjusting so that the magnetic 
field tends to be perpendicular to the axis of the filament. 




Fig. 5. — Number of stars (top), total stellar mass M* (mid- 
dle), and mass of the most massive star Mp (bottom) for all three 
runs as a function of free-fall time. In this figure and throughout 
the rest of this paper, we only count a sink particle as a star if 
it has passed the minimum merger mass of 0.05 M©, ensuring its 
permanence as the simulation proceeds. 



3.3. Fragmentation and Star Formation 

The most dramatic difference between the three runs 
is in the fragmentation. In all three cases, there is a pri- 
mary with a mass of about 23M0. In run BR, there is 
also a secondary star with less than IM© of material. In 
runs BI and HR, however, the filaments that feed the pri- 
mary object have fragmented into dozens of stars by the 
end of 0.6 free-fall times, with typical masses of O.2M0 
but ranging up to ~ 11M0. This filament fragmenta- 
tion takes place beginning around 0.2 to 0.3tff . By O.Stff , 
these stars have fallen into the central region and un- 
dergone significant N-body interactions with each other. 
After that time, the positions of the sinks in Figures 1 
and 2 no longer correspond to the places they were born 
- many of the sinks have been ejected towards the outer 
regions of the core. 

We summarize the properties of star particles in all 
three runs in Figures [5] and [6) Note that we only count 
a sink particle as a star once it has passed the mini- 
mum merger threshold of O.OSM©. Thus, the extra stars 
in runs BI and HR are not temporary objects that will 
eventually accrete onto the primary. While the exact 
value of this threshold is somewhat arbitrary, we point 
out that in runs BI and HR, there are a few dozen small 
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Fig. 6. — Fraction of total stellar mass that is in stars with mass 
less than m for all three runs at tff = 0.6. 
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sink particles that do not meet this threshold by 0.6 free- 
fall times, while in run BR there are none. Thus, we do 
not believe that the basic conclusion that fragmentation 
is dramatically suppressed in run BR compared to the 
others is sensitive to the exact numerical value of value 
of the minimum merger mass. 

One possible explanation for the difference between 
runs BR and HR is that extra fragmentation in run HR 
is due to disk fragmentation that is not present in the 
other runs because the magnetic field has removed much 
of the angular momentum from the central region. How- 
ever, this is not the case. From Figure 2, we can see that 
run HR has already undergone significant fragmentation 
in filaments well before the disk has grown large enough 
to fragment. In fact, most of the stars in run HR form at 
distances of a few thousand AU or greater from the cen- 
tral star - well outside the disk. Whatever the cause for 
the difference in fragmentation between the runs with a 
magnetic field and run HR, it is not due to the presence 
of a disk in one and not in the others. 

As discussed in section 2.2, although much of the frag- 
mentation in run BI is numerical in that it comes from 
filaments that collapse down to Pmax, it is still possible to 
choose Axl such that the fragment masses are roughly 
correct. From Equation (14), Pmax in run BI ranges from 
- 10-^"^ (for /3 ^ oo) to - 10-^^ g cm-3 (for (3 = 0.01), 
and so the density at which sink creation occurs in our 
simulations roughly mimics the density at which molec- 
ular gas can no longer cool efficiently. Thus, we expect 
that the fragmentation in run BI is qualitatively simi- 
lar to what would happen in massive cores if there was 
no protostellar feedback: the filaments would fragment a 
bit after they reached densities of ~ 10"-*^^ g cm~^, and 
one would end up with many more fragments than would 
be formed in the presence of radiative heating. Further- 
more, some of the protostellar properties in run BI do 
indeed appear to be converged. If we compare both the 
total mass in stars and the mass of the primary in run BI 
to the high-resolution version of BI at 0.2 tff , we find that 
they differ by only 6% and 5% respectively, over a fac- 
tor of 8 difference in resolution. Thus, while the number 
and mass distribution of the fragments in run BI are not 
converged, quantities that depend mainly on the overall 
accretion rate do seem to be. 

In addition to the fragmentation, we also find that the 
magnetic field slows down the overall rate of star forma- 
tion by abou t a fact or of about 3, consistent w i th Padoan,, 



fc Nordlund] ( |201l| ) and [Federrath fc Klessen| ( |2012l) At 
0.6 free- tall times, the total mass in stars in run BR is 
about ~ 2OM0, almost all of which is in the primary, 
compared to over 6OM0 in the run HR. Almost all of the 
"extra" star formation in the run HR has gone into stars 
other than the primary, which contains only ~ 40 % of 
the total stellar mass at 0.6 free-fall times. The mass 
of the most massive star, on the other hand, is approxi- 
mately the same in all three runs, probably because our 
initial conditions place the same amount of mass in po- 
sition to quickly collapse towards the center. Beginning 
at around 0.5tff, there is an increase in the rate of star 
formation in run HR as compared to the others. This in- 
crease is associated with the fragmentation of a filament 
formed in the outer region of the core that by ~ O.Stff has 
begun to form stars, as shown in the bottom panels of 



Figure 1. The relative timescales here are roughly what 
one would expect from inside-out collapse given our ini- 
tial conditions: for a power-law density profile with slope 
— 1.5, the ratio of the free-fall time at 0.75Rc to that at 
0.25i?c is about a factor of 1.5, which is approximately 
the delay we see here. Note that filament fragmentation 
in the outer regions of the core does not happen in either 
of the runs with a magnetic field - there, star formation 
only occurs close to the core center. We will discuss this 
difference further in section 4. 

We mention here as a caveat that our 10 AU resolution 
means that we cannot resolve any binaries closer than 
~ 40 AU, the accretion radius on one sink particle. Thus, 
we cannot rule out the possibility that the massive star 
present in run BR would in fact be massive binary with 
a separation of < 40 AU if we had higher resolution. 
However, even if that were the case, the fragmentation 
would still qualitatively different than in runs HR and 
BI, where we form dozens of stars with a masses that 
sample the fuh IMF. 

3.4. Thermal Structure 

The difference in fragmentation between runs BI and 
BR is expected, since it is well-established that radia- 
tive feedback in massive cores reduces fragmentation by 
raising the thermal Jeans mass of the collapsing gas (e.g. 



Krumholz et al. 2007a, 2010). The difference in frag- 



mentation between runs BK and HR, however, is more 
interesting. One possibility is that protostellar heating 
is somehow more efficient in the presence of magnetic 
fields. Figure [7] shows maps of the average temperature 
through a 5000 AU cube centered at the most massive 
star in runs BR and HR. We find that, contrary to this 
hypothesis, the heating in run HR is either similar to or 
slightly more widespread than in run BR, because accre- 
tion rates are higher in the absence of the field. This is 
not a dramatic effect, however. The total protostellar lu- 
minosity in run BR is typically smaller than that of run 
HR by only a factor of ^ 0.7. The temperatures, which 
in the optically thin limit scale like L^-^^, would be lower 
by only a factor of ~ 0.9. At 0.25tff, when the first frag- 
mentation in run HR occurs, the mean Tg in the 5000 AU 
cube around the primary is 63.3 K in run HR and only 
53.7 K in run BR, but despite the higher temperatures 
the gas in HR fragments while the gas in BR does not. 
So, the difference in the effectiveness of radiative heating 
between runs BR and HR cannot be responsible for the 
difference in fragmentation - it is too small and in the 
wrong direction. 

The difference, then, must be due to the direct support 
provided by the magnetic field in run BR. To quantify 
this effect, we define an effective temperature Teff by 



8^' 



(22) 



where n is the number of particles per unit volume. Ex- 
pressed in terms of /3, we find 



teff 



2/3 



(23) 



In other words, Teff is the temperature defined in terms 
of the thermal plus magnetic energy densities instead of 
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Fig. 7. — Maps of the average gas temperature, taken at the 
same times as Figures 1 and 2. The averages were taken along the 
X direction through a 5000 AU cube around the most massive star. 
Run BR is on the left and run HR is on the right. 



just the thermal energy density. Note that this is actu- 
ahy more closely related to our criterion for creating a 
sink particle than the gas temperature because we have 
included the magnetic energy in defining the magnetic 
Jeans number (see the Appendix). While the concept 
of an effective temperature is clearly an oversimplifica- 
tion - for one, the magnetic field does not resist collapse 
isotropically the way thermal pressure does - we find that 
it is helpful in understanding our simulation results. 

In Figure |8j we summarize the combined temperature 
and magnetic field structure of the cores in runs BR, BI, 
and HR. In the two right panels, we plot the total mass 
in each p — Tg bin for runs BR and HR over a series of 
time snapshots. In the two left panels, we instead use p — 
Teff bins for the two runs with magnetic fields. The top 
row of the figure merely summarizes our initial condition. 
Although the core temperature starts at precisely 20 K 
in all three runs, the cylindrically symmetrical magnetic 
field profile means that there are a range of magnetic 
field strengths, and thus Teff covers a range of values. 
The blue diagonal lines represent the threshold at which 
the code lays down a sink particle. Thus, there can be no 
gas in any of the runs to the right of this line - any cell 
that exceeds this threshold has some of its gas converted 
into sink particles until it no longer violates the MHD 
Truelove criterion. This line is suppressed in the third 
column, because in the presence of magnetic fields, there 
is no single density at which sinks are created for a given 
temperature (see Equation 14). 

We can get a sense of whetner star formation is taking 
place from these plots by looking at whether there is any 
gas close to crossing this threshold. In run BR, there 
is hardly any gas close to the densities required for sink 
formation. Runs BI and HR, on the other hand, have sig- 
nificant amounts of gas close to that threshold by around 
0.2 to O.Stff . The phase diagram for run HR, in particu- 
lar, bears a number of "finger" features that correspond 
to gas that is all at one Tg , but that stretches over a range 
of densities approaching that required for sink formation. 
These features are most prominent at O.Stff, but are vis- 
ible before and after as well. The is precisely the time at 
which the main filament in run HR has broken up into 
a number of gravitationally unstable "beads" , which col- 
lapse down until they form sink particles. The "fingers," 
then, correspond to gas in these beads that is collaps- 
ing isothermally, albeit at higher temperatures than the 
initial 20 K, with the precise value determined by the 
distance from the bead to the central protostar. This 
collapse is isothermal because the temperature changes 
on the evolution timescale of the most massive protostar, 
which for our problem is tff, while the timescale for lo- 
cal gravitational collapse in the bead is must faster. In 
contrast, we do not see this behavior in run BR, because 
a combination of magnetic and thermal support has ren- 
dered the main filament in that run stable against grav- 
itational collapse at a density much higher than the sink 
creation value. 

In one sense. Figure [S] restates what we already know 
- there is much fragmentation in runs BI and HR and 
hardly any in run BR. However, this plot can also help 
us untangle the effect of the magnetic field and the radia- 
tive feedback by telling us in which regimes each effect is 
more important. By comparing runs BR and BI, for in- 
stance, we can see that the primary effect of the radiation 
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is to heat up the relatively dense regions in the core - i.e. 
to move material greater than about lxlO~^^ g cm~^ up 
in the plot and away from the sink formation threshold. 
Alternatively, the slope of the (p — Teff ) phase diagrams 
for runs BR and BI show that the magnetic field is most 
effective at raising Teff at low density. Hence, we can 
begin understand that the reason the combination of the 
B field and the radiative feedback is more effective at 
suppressing fragmentation than either considered in iso- 
lation is that they are effective in different regions, with 
the magnetic field mostly helping to support (or, at least, 
to slow the collapse of) material in the diffuse, outer parts 
of the core, and with radiation most effective in the dense 
material that is close to the central protostar. 

4. DISCUSSION 

4.1. Why do Magnetic Fields and Radiation Suppress 
Fragmentation? 

We would like to understand the suppression of frag- 
mentation in run BR in terms of the mass-to-flux ra- 
tio The average /j.^ for the entire core is 2, but, 
because the core is centrally concentrated, it is greater 
through flux tubes passing near the center and lower 
through flux tubes passing through the diffuse, outer re- 
gions. It is illustrative to do the following analysis on 
our initial conditions: take the initial spherical region 
and exclude a cylindrical region of radius Rz concentric 
with the sphere and extending through the entire do- 
main. Then, compute /i^,^^ the mass-to- flux ratio in the 
remaining region, ja^^z is 2 when R^ = and mono- 
tonically drops to when Rz = Rc How quickly jii^^z 
drops off with Rz will give us a rough estimate of where 
we can expect the core to be subject to fragmentation. 
We find that, by a radius of Rz ~ 0.73 Rc^ /i^,^ has 
dropped below 1, meaning that the region external to 
that cylindrical radius (corresponding to approximately 
32% of the core volume and 19% of the mass) should be 
fairly well-supported against collapse. Furthermore, the 
point at which /j^^^z has dropped to 1.5 is at only 0.44 
Rc, meaning that '^72% of the core volume and ~53% 
has a mass-to- flux ratio below that value. While struc- 
tures with a mass-to- flux ratio of 1.5 are supercritical and 
should collapse, they will still collapse more slowly than 
in the absence of the magnetic field, giving the radiative 
feedback more time act. This effect is not dramatic; the 
effect of the magnetic pressure force in the virial theo- 
rem is to dilute g:ravity al on^ the field lines by a factor of 
(1 — /ij^) ( Shu fc Li|l997 ), so that structures that are su- 
percritical Dy~aTac^r~oi 1.5 collapse approximately half 
as quickly as structures with of infinity, and even the 
core as a whole collapses about 75% as fast at /i<|, = 2. 

Thus, even inside a supercritical core, the magnetic 
field can slow collapse in specific sub-regions with mass- 
to- flux ratios below unity, or even halt it altogether, and 
this mechanism is more effective at suppressing fragmen- 
tation in the outer regions of the core. We can see this 
effect operating in the bottom row of Figure 1 - while 
most of the star formation in run HR takes places to- 
wards the central region of the core, by 0.6 free-fall times 
we have begun to see signs of fragmentation of a filament 
in the outer regions as well. This does not take place in 
either of the runs with magnetic fields, including the one 
with no radiation, so this cannot be a radiative effect. 



Rather, it is due to the ability of the magnetic field to 
effectively suppress fragmentation in the outer regions of 
the core. To quantify this, we plot in Figure [9] the dis- 
tribution of r*, which measures how far away each star 
was from the central massive star at the time it formed. 
We compute this quantity for every star (other than the 
first) that forms over the history of each simulation. In 
run BI, most of the fragmentation takes place at distance 
of ~ 2, 000 AU, and there are no stars that form at a dis- 
tance greater than ~ 5, 000 AU from the central object. 
In run HR, however, the most likely value of r is roughly 
5000 AU, while a significant fraction (about 1/3) of the 
stars form at distances of 10,000 AU or greater. Star 
formation at such large radii is completely suppressed by 
the magnetic field, even without heating effects. 

The following picture thus emerges: magnetic fields 
work to suppress fragmentation in the outer regions of 
the centrally-concentrated cores, either by slowing it 
down or halting it altogether. If they halt it altogether, 
then fragmentation is confined to the central region, 
where radiative heating is most effective. If magnetic 
fields merely slow fragmentation at large radii, then they 
still allow radiative heating more time to "win" by heat- 
ing up filaments to the point at which they are too warm 
to collapse further. The combined result is that the mag- 
netic field and the radiation are together far more effec- 
tive at suppressing fragmentation than either process in 
isolation. Our work suggests that typical massive cores, 
which are centrally concentrated and have /i*^, ~ 2, do not 
fragment strongly, as one would expect from the corre- 
spondence between the core and initial stellar mass func- 
tions. 

Finally, while our simulations are based on ideal MHD, 
we do not expect that non-ideal effects will dramatically 
alter our conclusions. Ohmic dissipa tion, ambipolar dif- 
fusion, and reconne ction diffusion ( [Santos- Lima et al. 
2010 Lazarian|2011 ) are all capable of increasing /i<|, , but 
they also all most important at high densities and/or re- 
gions where the magnetic field lines are most bent. Those 
are precisely the regions where the magnetic field is least 
important for the suppression of fragmentation, because 
they are all primarily associated with the inner region of 
the core where radiative heating is most effective. 

4.2. Comparison to Commergon et. al. 

"Commercon et al. (2011) have also performed a set of 
radiation-magnetohyorbdynamic simulations of massive 
core collapse using similar initial conditions to the ones 
considered here. They also found a synergistic effect be- 
tween the radiative heating and the magnetic field, where 
the two effects in tandem lead to much less fragmentation 
than either considered in isolation. However, their work 
differs from our own in a few key respects. Most signifi- 
cantly, they have focused on the early stages of collapse, 
up to just past the point at which the first hydrostatic 
core forms, while we have focused on what happens to 
the remainder of the core after the first protostar has 
undergone second collapse. Thus, while we have both 
identified mechanisms by which a combination of radia- 
tive and magnetic effects suppress fragmentation, these 
mechanisms have different underlying causes and mani- 
fest themselves at very different times. 

This difference in emphasis stems from our different 
recipes for representing protostellar feedback. In this 




Fig. 8. — Left two panels - phase diagrams showing the amount of mass in each p — T^q bin at different times, for runs BR and BI. Right 
two panels - the same, but with p — Tg bins for runs BR and HR. The snapshots are taken at the same times as the above figures. The 
islands of low density material at Tg ~ 10^ K and Teff ~ IC^ K correspond to gas in the ambient medium and should be ignored. 
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Fig. 9. — Histograms of star formation distance r* for runs BI, 
HR, and BR. Here, r* is the distance each star was from the most 
massive star when it formed, computed for every star that forms 
over the entire history of each simulation. In run BR, there is only 
one secondary fragment, which forms a distance of r* ~ 3600 AU 
from the primary. 

work, we use sink particles to represent material that 
has collapsed to densities higher than we can follow on 
the grid, and compute a luminos ity Lj for each p article 
according to a sub-grid model. |Comme rcon et al. ("2011), 
on the other hand, do not use sink particles, but instead 
employ higher resolution (2.16 AU), which allows them 
to follow the formation of a hydrostatic core. The feed- 
back from the accretion shock on to this core can then be 
computed on the grid. They find the radius at which this 
shock releases its energy depends on the core's magnetic 
field, since strong magnetic braking allows smaller first 
cores to form. 



The approach of Commercon et al. (2011) has the ad- 
vantage of self-consistency, l^urtherniore, it is probably 
more accurate than our technique at representing radi- 
ation from particles before the second core has formed. 
In fact, we do not include any radiative feedback until 
the sink particle mass exceeds O.OIM©, and second col- 
lapse generally takes place at a few times that value, so 
except for a very short period of time, we ignore this radi- 
ation altogether. Also, because their resolution is higher 
and they do not have to handle sink particle mergers, 
they can resolve binaries that we cannot. However, their 
method has the downside that it cannot model the effect 
of the much larger (by a factor of ~ 100) accretion lumi- 
nosities that occur after second collapse. Additionally, 
the lack of sink particles s everely limits the integration 
time for Commergon et al. (12011), since they cannot fol- 
low collapse past the point wEere they fail to resolve the 
Jeans density. This limited them to running for only a 
few percent of a free-fall time after the first hydrostatic 
object formed. In contrast, with our initial conditions 
the first core forms almost immediately, and we find that 
most fragmentation does not occur until around 20% to 
30% percent of a free-fall time past that point. In a 



sense, our simulations pick up where those of |Commercon| 
et al. (^2011) left off, in that our simulations begin with a 
itrair 



Sub-grid luminosity models have problems of thei r own 
related to unresolved binarity, as pointed out by [Bate 
([2012j). As mentioned in section 3.3, our use of sink par- 
ticles means that we cannot resolve any binaries closer 
than 40 AU, and we cannot rule out the possibility that 
the central massive star in our simulations in fact rep- 
resents an unresolved binary. However, for accretion 
luminosity-dominated stars, it matters little whether a 
sink particle represents a single star or a binary too 
tight to be resolved, because the energy released per unit 
mass accreted onto low-mass protostars is nearly inde- 
pendent of the stars' masses (jKrumholz 201 1|. On the 
other hand, stars' internal luminosity scales with mass 
as roughly M^-^, meaning that in the worst case where 
a sink particle should in fact represent an equal mass 
binary, the internal luminosity is overestimated by a fac- 
tor of 2^-^ = 5.7. While this is a potential concern, the 
internal luminosity does not become comparable to the 
accretion luminosity in our calculations until about 0.3tff , 
and by that time there are already clear differences in the 
fragmentation between runs HR and BR. Moreover, the 
alternative of not including a sub-grid luminosity model, 
as in Bate (2012), is far worse. Without such a model one 
omits both the accretion luminosity onto the stellar sur- 
face and the larger internal luminosity, and the resulting 
error is many orders of magnitude. 

Finally, we mention one last diff erence between our 
work and |Commercon et al. ( 2011[ ): their initial condi- 
tions contained much less kinetic energy than our own, 
with avir = 0.2 versus avir = 2.3. This could explain 
why we see a sma l l disk in our high resolution run, and 
Commercon et aP (2011 1 do not. If the ir /3rot ^ 0.02avir5 
as implied by|Burkert ^ Bodenheimer] (^000 ) , then they 
would have /3rot ^ 0.004, smaller th an our own by a fac- 
tor of 3. On the other hand, Seifried et al. (2012) had 
P = 0.04, higher than ours by a factor of 4 again. I'hus, 
the sequence of disk sizes seen in ou r papers, ranging 

40 AU (us 



2012D to 



from - 100 A U (iSeifried et al. [ 

to unresolved ("Commercon et a i.||2011 ), could simply be 
a consequence of different amounts oi angular momen- 
tum in the co r es. I t is possible, however, that if Com-] 
mercon et al. ( |2011^ ) extended their simulation to later 
times, that tney too would begin to resolve a disk in 
their ja^ = 2 run. We conjecture that this disk would be 
smaller than ~ 40 AU in radius. 

4.3. Where is Fragmentation Suppressed? 

An interesting question is: in what range of the S — /i<^ 
pa rameter space is fragm entation weak? As discussed 
inlCrutcher et al. (2010), although the average molecu- 



centrally concentrated core with one protostar that very 
quickly undergoes second collapse. 



lar cloud core is marginally magnetically supercritical, it 
by no means follows that there are no cores with weak 
magnetic fields. The Bayesian analysis presented in that 
paper suggests the distribution of field strengths is quite 
fiat, such that there may be many cores where the field is 
significantly weaker than the ones discussed here. In fact, 
for cores like ours with a mean density of nn = 2.4 • 10^ 
cm~^, their result suggests that the total magnetic field 
strength should be evenly distributed between ^0.0 mC 
and ~ 3.4 mG, roughly twice the value considered here. 
This implies that about 25% of cores like the ones in this 
paper would have values of /i«^> of 4 or greater. Our work 
suggests that there should be a tendency towards greater 
fragmentation in massive cores with such weak magnetic 
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fields, with such cores being more likely to form clusters 
rather than isolated massive stars or binaries. A recent 



set of millimeter observations (Palau et al. 2013) stud- 
ied the fragmentation of 18 massive cores with < 1000 
AU resolution, and found that ~ 30% showed no signs of 
fragmentation, while 50% did. They propose that vari- 
ation in the magnetic field strength may be responsible 
for the determining the fragmentation, but confirmation 
of this view will have to wait for follow-up observations 
of the field. 

Furthermore, the recent observations oflButler & Tanl 



(2012) found a typical massive core surface density of 
~ 0.1 g cm~^, over a factor of 10 lower than the 2 g 
cm~^ cores considered here. These cores are below the 
surface density thres hold for massive star formation ~ 1 
g cm~^ identified in Krumholz & McKeel (2008), which 
ignored magnetic fields. Could magnetic fields play a role 
in lowering the threshold for massive star formation? We 
plan to address these questions in future work. 



5. CONCLUSIONS 

We have presented a set of 3D, R-MHD simulations of 
the collapse of isolated, magnetized, massive molecular 
cloud cores that attempt to isolate the effects of the mag- 
netic field and of the radiative feedback on core fragmen- 
tation. We find that the magnetic field and protostellar 
radiation can combine to largely suppress fragmentation 
throughout the core, so that the simulation that includes 
both magnetic fields and radiation results in only a sin- 



gle binary star system, while the runs that exclude either 
effect are subject to far more fragmentation. The expla- 
nation for this behavior is that magnetic fields and radia- 
tive heating are effective in different regimes, so that that 
each effect influences gas that the other misses. We find 
that massive cores with typical magnetic field strengths 
likely collapse to form single star systems, as suggested 
by the observed relationship between the CMF and IMF. 
We have also reproduced the result found by other re- 
searchers that Keplerian disks can form in the presence 
of magnetic fields with ~ 2, provided that turbulence, 
which results in a misalignment between the magnetic 
field and angular momentum vectors, is present. 
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APPENDIX 
A. THE MHD TRUELOVE CONDITION 
Simulations of isothermal, self-gravitating systems are subject to artificial fragmentation unless the Jeans length, 

.2\ 1/2 



Aj = 



'Gp 



is resolved by a sufficiently large number of grid cells. 



Ax 



th 



th, max 7 



(Al) 



(A2) 



where Ax is the width of a grid cell fTru elove et al.||l997 ). We have added the subscript "th" to the Jeans number 
J to indicate that it is for the purely thermal case, m which there is no magnetic field. For the case they studied, 
[TVuelove et al.| ( |1997| ) found that Jth, max — 0.25 was adequate to suppress artificial fragmentation, but in general it is 
a problem-dependent quantity. For AMR simulations, the Truelove condition is one of the criteria used to increase the 
refinement; for simulations in which sink particles are used as a sub-grid model for protostars, the Truelove condition 
is often used to determine when sink particles should be introduced — i.e., whenever the density exceeds 



Pn 



^'"^th,max^8 

GAx2 



(A3) 



Magnetic fields suppress fragm entation and therefore should allow one to defer refinement or the introduction of 
sink particles to higher densities. [Federrath et al. (2010) included the effects of the magnetic field in their refinement 
criteria, but it was used as a supplement to the thermal Truelove criterion, not as a replacement. For introducing sink 
particles, they did require that the total energy of a control volume be negative. 

To generalize the Truelove condition to include magnetic fields, we b egin w ith the expression for the maximuni mas s 
of an isothermal, magnetized cloud derived by |Mouschovias fc Spitzer (1976), as generalized by Tomisaka et al. (1988), 



Mcr = I.ISMbe 



-3/2 



(A4) 
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where 



is the Bonnor-Ebert mass and 



Mbe = 1.18 



(A5) 



(A6) 



is the magnetic critical mass. (Tomisaka et al. (|1988 ) found that a factor 0.17 fit their numerical results better than 
l/(27r), but we adopt the latter lor simplicity.) Using the alternative form for the magnetic critical mass, M^, which 
is defined by 



Equation (A4) can be expressed as 



3/2 



(A8) 



dBertoldi fc McKee| ( |l99"2| , who wrote Mj for the Bonnor-Eber t mass). Evaluation of Mb gives Mb /Mbe = 0.76/3 



-3/2 



We define the critical radius by ( Mouschovias fc Spitzer|[l976 ) 
and then obtain 

Rcr = 0.39Aj ( 1 



4 



0.74 



1/2 



(A9) 



(AlO) 



We denote the critical radius in the absence of a magnetic field {{3 oo) as i^cr,th- As expected, we see that 

2i?cr,th — Aj. 

In the non-magnetic case, the Jeans number is Jth = Ax/Xj ex Ax/i?cr,th- We then generalize the Jeans number to 
the MHD case by writing 

f = (All) 

^th ^cr 



SO that 



T,. 0.74 V^' 

Jth = J ( 1 + 



(A12) 



Sinc e the same relation applies to the maximum Jeans numbers, the MHD Truelove condition follows from Equation 
(A3|: 

Pmax = „9 ( 1 H — ) • (A13) 



We note that if one expresses the Jeans length in terms of the energy density, u = |pc^, as Aj = (27rii/3Gp^)^/^ and 
then adds the magnetic energy density into one obtains the same result as in Equation (A13) except that the factor 



0.74 is replaced by |, a negligible difference. Our result is thus very similar to the approach advocated by Federrath| 



et al. (2010). The advantage of the present derivation is that it is directly tied to the maximum stable mass. 

Magnetic fields c an halt collapse perpe ndicular to the field, but they have no effect on gravitational instability 
parallel to the field ( Chandrasekhgir|[l961 ) . Gas that collapses along the field lines has a thickness 



where S is the surface density, so that 



H 



H 



v/2 



Po 



Aj, 



0.45 



Ax. 



(A14) 



(A15) 



In fact, a self-gravitating sheet cannot become thinner than 2 Ax, since a single layer of cells cannot exert a vertical 
gravitational force inside the layer. Hence, if it is important to follow the internal dv namic s of gravitationally stable 
sheets, one should maintain 2 Ax < corresponding to Jth ^ 0.25 from Equation (A15). This criterion does not 
apply to gravitationally unstable sheets (J > Jmax), since they will either be refined or replaced by sink particles. 
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Fig. 10. — Logarithm of the column density (normahzed to p(0) -^o) through the filament in each run at the point at which the maximum 
density has reached 2 x 10^ times the initial value. Top Left - Jmax = 0.125. Top Right - Jmax = 0.25. Bottom Left - Jmax = 0.375. 
Bottom Right - Jmax = 0.5. The top panels have only one filament, while the bottom two show clear signs of artificial fragmentation. Each 
image shows a region of size ~ 0.07 Iq. 



To test the MHD Truelove criterion, we carry out the same test used by Truelove et al. (1997|), but with the addition 
of an initiahy uniform magnetic field. We begin with a cubic, periodic box ot size Iq tilli 
has a spherically symmetric, Gaussian density profile 



ed with isothermal gas that 



p(r) = p(0) exp 



(A16) 



Here, p(0) is the central density and ri is a characteristic fall-off radius, which have taken to be 0.48 Iq. To this 
background density we also added a m = 2 azimuthal density perturbation with an amplitude of 10%. We report our 
simulation results in units that have been normalized by p(0), io^ and the central-density free-fall time, given by 



32Gp(0) ' 



(A17) 



We have also set the core in initial rotation about the z axis with an angular velocity uj such that ujtf[ = 0.5 3. This 



value was chosen so that the core would make approximately 1 rotation in 12 free-fall times, as in the original [Boss 
[Bodenheimer" ("IDTQl) version of this problem. Finally, we have imposed an initially uniform magnetic field field pointing 
m the z direction, with a magnitude such that M/M^ ^ 2.25. This field is strong enough to alter the morphology of 
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the collapse, but not strong enough to halt it completely. This value of the mass-to-flux ratio is also comparable to 
that observed in real star-forming regions. The initial plasma beta /3o is ~ 1.57 at the center of the domain and as 
low as 0.19 at the edges. 

As in the no n- magnetic case, the gas first collapses into a sheet and then into a filament, where the field is normal to 
the filament. The plasma beta in the sheet before it begins to collapse is of order unity. Since the magnetic energy and 
the gravitational energy are both independent o f the radius of the filam ent, an isothermal filament with a supercritical 
mass-to-flux ratio will col lapse indefinit ely (In utsuka fc Tsuribe||2QQl ). To see this, note that the gravitational energy 
per unit length is —Gmj ( )Fiege fc PudritZ||2QQQ ), where is the mass per unit length, whereas the magnetic energy 
per unit length is of order 7rr^5^ oc where is the magnetic flux per unit length. (The exact expression for the 
magnetic energy depends on the structure of the field inside the filament.) Thus, the ratio of the two forces depends 
only on the mass-to- flux ratio per unit length, which is constant in ideal MHD. Note, however, that the scaling of the 
magnetic field is only valid if the field is perpendicular to the filament axis. 

We have run four versions of this problem, each time using the AMR capabilities of our code to impose a different 
^max- The results are shown in Figure [lOj Because the free-fall time is a function of density, more or less well-resolved 
simulations of this problem will not be at the same stage of development at the same simulation time. We instead 
compare the runs at the point where they have all reached approximately the same maximum density of 2 x 10^ p(0). 
The simulation times at which the peak density reaches this value range from about 3.79 in the best-resolved case 
to 4.02 tff in the worst. We find that, as in the pure hydro version of this problem, Jmax = 1/4 is sufficient to halt the 
onset of artificial fragmentation. The maximum thermal Jeans number, on the other hand, is larger then Jmax by a 
factor of ~ 2 in these runs. Refining on the less stringent magnetic Jeans number seems sufficient to accurately follow 
the collapse of magnetized gas for this problem, although other physical processes, such as B-field amplification via 
dynamo action, may require a higher resolution ( Federrath et al.||2QlT| . 
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